A first stab at a diagnostic surface temperature implementation #278

simone-silvestri · 2024-11-26T10:33:57Z

This PR implements two types:

BulkTemperature
SkinTemperature

The first type is supposed to represent the default option we have know: the surface temperature is provided and equal to the temperature of "surface" model (either sea ice or ocean). This temperature is fixed and not iterated over.

The second type represents a diagnostic temperature that is diagnosed by a flux balance between the

top turbulent surfaces (dependent on the surface temperature and iterated by the solver)
top radiative fluxes (prescribed)
outgoing longwave (calculated from the surface temperature in the previous iteration)
internal fluxes: provided from the SkinTemperature, for the moment only a simple diffusive flux with prescribed diffusivity leading to an explicit expression for temperature, but they can become more complicated.

In this second case the temperature is solved for in the flux_balance_temperature(::SkinTemperature, T_ocean, net_external_flux) function

These types are provided to the similarity_theory as a new field named surface_temperature_type. The surface temperature is then stored in a T_surface field which is added to similarity_theory.fields

There is a bit of a pitfall in this implementation: it is difficult to implement temperature-dependent radiative properties, because radiative properties (that depend on i and j and the grid) are passed already computed to the solver.
I would be against passing i, j and the grid to the surface flux solver (it would entangle the solver with oceananigans so it would be difficult to swap it out with a SurfaceFluxes.jl solver), so we could rethink albedo and emissivity computation.

…into ss/one-degree-take3

simone-silvestri · 2024-11-26T14:59:55Z

A very crude first calculation in experiments/single_column_experiments/os_papa_surface_temperature.jl shows the difference between prescribed and diagnostic surface temperature. In the latter case, the internal fluxes are just a very simple diffusive flux with

$$\kappa \frac{T_s - T_o}{\Delta z}$$

where $\kappa = 10^{-2}$, $T_s$ is the incognita, $T_o$ is the temperature of the first ocean cell, and $\Delta z$ is half the upper cell spacing.
The solid lines are the results of DiagnosticSurfaceTemperature and the dashed ones are the respective variables calculated using PrescribedSurfaceTemperature

single_column_profiles.mp4

There is a large difference in the ocean surface temperature (about 1 degree K), but very little in the fluxes.

simone-silvestri · 2024-11-26T15:12:36Z

However, the surface diffusivity predicted by CATKE in the above video is $O(1)$, much larger than the $\kappa = 10^{-2} m^2s^{-1}$ that I prescribed. If I use $\kappa = 0.5 m^2s^{-1}$ there is even less difference in the surface temperature and barely any in the heat fluxes.

single_column_profiles.mp4

I need to test this on the whole surface to see if there are maybe hot or cold spots that show more difference. As it stands, the net difference in heat flux is less than 1 $Wm^{-2}K^{-1}$. Difficult to assess the impact on a 10-day timescale.

Probably for sea ice where timescales are faster this might matter more.

simone-silvestri · 2024-11-26T15:46:49Z

There is actually quite a large difference in net heat flux: This is the same computation as in the generate_surface_fluxes.jl but done also with a DiagnosticSurfaceTemperature with a diffusivity of 0.1.
Here I am showing the difference of net heat flux (Q_prescribed - Q_diagnostic) on the left and surface temperature (Ts_prescribed - Ts_diagnostic) on the right.
The surface temperature difference is upward of 2 degrees and the net heat flux has differences upward of 40 Wm2K-1 (which is quite large). In practice, a diagnostic temperature seems to have much more extreme fluxes, i.e heating regions (where solar radiation peaks) are heating much more and cooling regions are cooling much more. There might be a bug in the solver though, because some regions NaN.

simone-silvestri · 2024-11-26T16:19:58Z

There were some gross bugs in the implementation (including the sign of the fluxes). Correcting the implementation the results look like this (same plot as before)

So in general a diagnostic temperature seems to be mitigating (slightly) the surface fluxes as it is intuitive it should do.
Surprisingly enough, it seems that the solver converges in average more quickly with a Diagnostic surface temperature

src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl

codecov · 2024-12-02T16:19:10Z

Codecov Report

Attention: Patch coverage is 0% with 67 lines in your changes missing coverage. Please review.

Project coverage is 0.00%. Comparing base (0fff4fa) to head (08be1bb).
Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
...sRealmFluxes/similarity_theory_turbulent_fluxes.jl	0.00%	32 Missing ⚠️
...aIceModels/CrossRealmFluxes/surface_temperature.jl	0.00%	20 Missing ⚠️
...Models/CrossRealmFluxes/atmosphere_ocean_fluxes.jl	0.00%	9 Missing ⚠️
...rc/OceanSeaIceModels/CrossRealmFluxes/radiation.jl	0.00%	6 Missing ⚠️

Additional details and impacted files

@@          Coverage Diff          @@
##            main    #278   +/-   ##
=====================================
  Coverage   0.00%   0.00%           
=====================================
  Files         33      34    +1     
  Lines       1876    1910   +34     
=====================================
- Misses      1876    1910   +34

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl

Co-authored-by: Gregory L. Wagner <[email protected]>

simone-silvestri · 2024-12-05T07:07:40Z

Shall we merge?

glwagner · 2024-12-05T14:16:32Z

There was a question about recomputing the saturation specific humidity. Was this changed and did it change the posted results?

glwagner · 2024-12-05T14:17:26Z

It's a very nice result that it converges faster by the way.

simone-silvestri · 2024-12-05T14:21:06Z

There was a question about recomputing the saturation specific humidity. Was this changed and did it change the posted results?

Yep I changed it and indeed, the difference between the different solvers is much larger now in the order of 0.2 degrees instead of 0.02. I will post here another test showing the differences

glwagner · 2024-12-05T14:34:17Z

experiments/single_column_experiments/os_papa_surface_temperature.jl

I'd be ok to just have one Papa simulation, which uses the diagnostic surface temperautre (no need to keep the other one)

glwagner · 2024-12-05T14:35:42Z

src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl

+#
+# Where the LHS is the internal diffusive flux inside the ocean (within the boundary layer of thickness δ) 
+# plus the (semi-implicit) outgoing longwave radiation and the RHS are the remaining atmospheric and radiative fluxes
+# provided explicitly.


I appreciate the comment, but which method was chosen?

right, the fully explicit one. Should I remove the comment? Or maybe just update it

I think just add a remark about what the code does now; the linearization strategy is a potential optimization that isn't taken yet

or you can remove it :-D

glwagner · 2024-12-05T14:38:07Z

src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl

+                                               AtmosphericThermodynamics.Liquid())
+
+    𝒬ₛ = AtmosphericThermodynamics.PhaseEquil_pTq(ℂ, 𝒬₀.p, θ₀, q₀)
+    q₀ = AtmosphericThermodynamics.vapor_specific_humidity(ℂ, 𝒬ₛ)


You don't want to call this again, because the "seawater" part is important --- we account for the h2O mole fraction there. Also, we may be using a different saturation parameterization than AtmosphericThermodynamics (this is to comply with OMIP protocol)

hmm ok, if I remove that line the results are different (regression tests do not pass). Isn't ok that the seawater part called above?

weird, not sure I understand what is going on then

I mean, were the regression tests also incorrect?

Or maybe I am confused. Does vapor_specific_humidity(ℂ, 𝒬ₛ) compute a new specific humidity or merely extract it?

I think I'm confused sorry, I guess vapor_specific_humidity(ℂ, 𝒬ₛ) just extracts it. It says nothing about saturation there.

But I am weirded out that it changes the result. What is the difference between the two q0?

then probably yes, because in the main code we do:

ClimaOcean.jl/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl

Lines 380 to 382 in 0fff4fa

q₁ = AtmosphericThermodynamics.vapor_specific_humidity(ℂ, 𝒬₁)

q₀ = AtmosphericThermodynamics.vapor_specific_humidity(ℂ, 𝒬₀)

Δq = q₁ - q₀

where 𝒬₀ holds the previously calculated q₀

Ok, on main what we do is compute q0, and then we assemble it into a thermodynamic / dynamic state. Then we pass the whole state to the solver. Within the solver, I guess it uses vapor_specific_humidity to extract q0 --- or at least, that's what I thought. This result suggests though that somehow q0 is changed?

I think it is changed in the vapor_specific_humidity function. It is kind of tough to follow the thread, but it looks like the vapor_specific_humidity builds a PhasePartition type to calculate the humidity vapor fraction associated with the thermodynamic state

glwagner · 2024-12-05T14:56:25Z

There was a question about recomputing the saturation specific humidity. Was this changed and did it change the posted results?

Yep I changed it and indeed, the difference between the different solvers is much larger now in the order of 0.2 degrees instead of 0.02. I will post here another test showing the differences

Can you also show that the difference depends on the grid spacing at the top? It should right?

simone-silvestri · 2024-12-09T12:02:50Z

Can you also show that the difference depends on the grid spacing at the top? It should right?

It does not for the moment because the boundary layer thickness δ is a parameter passed to the SkinTemperature type.
In practice, the first cell is assumed to be within the well-mixed region, so T represents the (mixed) cell average, and δ is the boundary layer at the top of the mixed region. These are the differences with the default parameters for SkinTemperature (κ = 0.01 and δ = 1.0) (BulkTemperature - SkinTemperature)

These parameters probably require calibration.

glwagner · 2024-12-09T16:16:35Z

Can you also show that the difference depends on the grid spacing at the top? It should right?

It does not for the moment because the boundary layer thickness δ is a parameter passed to the SkinTemperature type. In practice, the first cell is assumed to be within the well-mixed region, so T represents the (mixed) cell average, and δ is the boundary layer at the top of the mixed region. These are the differences with the default parameters for SkinTemperature (κ = 0.01 and δ = 1.0) (BulkTemperature - SkinTemperature)

Hmm right, maybe it shouldn't depend on the grid spacing. That's helpful actually.

simone-silvestri and others added 14 commits November 22, 2024 14:15

this works best?

3114f74

One degree simulation take 4

9101d57

add a biharmonic viscosity

5833d19

Merge branch 'main' into ss/one-degree-take3

96a64a2

infer maximum delta t

0aa74dc

Merge branch 'ss/one-degree-take3' of github.com:CliMA/ClimaOcean.jl …

ea8708f

…into ss/one-degree-take3

Merge remote-tracking branch 'origin/main' into ss/one-degree-take3

3a77b67

changes

7c87f2b

this sohuld work

3e3b77d

using CUDA

d7ac6fa

tests NaN

9eaa668

test pass

6d094b1

regression test pass!

51328c8

regression tests pass!

33f69ab

simone-silvestri marked this pull request as ready for review November 26, 2024 12:50

simone-silvestri added 3 commits November 26, 2024 15:13

now it works

d1a3caf

reverted to false?

1f54525

the complete thing

3513848

chnage defaults

f754f7e

simone-silvestri changed the title ~~A first stab to diagnostic surface temperature~~ A first stab at a diagnostic surface temperature implementation Nov 26, 2024

switched a couple of variables

e6d8829

simone-silvestri added 4 commits November 26, 2024 17:29

remove comment

4e8faea

done better

af28e4e

better

9b32122

pass also gravity

d2f3fea

glwagner reviewed Nov 27, 2024

View reviewed changes

src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl Outdated Show resolved Hide resolved

simone-silvestri added 5 commits December 2, 2024 15:15

remove tyep instability

9313373

better syntax

542efa1

&& -> &

0030b9d

another bugfix

60ad704

Merge branch 'main' into ss/new-fluxes

13622fc

simone-silvestri requested a review from glwagner December 2, 2024 14:29

simone-silvestri added 4 commits December 2, 2024 15:33

do not change OceanSimulations.jl here

c8db739

added a test for SkinTemperature zero fluxes

8616feb

import types

54dd6c1

fix tests

fc8ddc0

glwagner reviewed Dec 2, 2024

View reviewed changes

src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl Outdated Show resolved Hide resolved

simone-silvestri and others added 2 commits December 2, 2024 18:26

Update src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl

2ea71e3

Co-authored-by: Gregory L. Wagner <[email protected]>

Merge branch 'main' into ss/new-fluxes

c225076

glwagner mentioned this pull request Dec 5, 2024

Refactor turbulent flux solver to have a DiagnosticSurfaceTemperature mode #162

Closed

glwagner reviewed Dec 5, 2024

View reviewed changes

simone-silvestri added 2 commits December 5, 2024 15:41

change comment + only one papa simulation

f915ed3

add comment

08be1bb

glwagner approved these changes Dec 9, 2024

View reviewed changes

glwagner merged commit b73c67c into main Dec 9, 2024
19 checks passed

simone-silvestri deleted the ss/new-fluxes branch December 9, 2024 17:47

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

A first stab at a diagnostic surface temperature implementation #278

A first stab at a diagnostic surface temperature implementation #278

simone-silvestri commented Nov 26, 2024 •

edited

Loading

simone-silvestri commented Nov 26, 2024 •

edited

Loading

simone-silvestri commented Nov 26, 2024 •

edited

Loading

simone-silvestri commented Nov 26, 2024

simone-silvestri commented Nov 26, 2024 •

edited

Loading

codecov bot commented Dec 2, 2024 •

edited

Loading

simone-silvestri commented Dec 5, 2024

glwagner commented Dec 5, 2024

glwagner commented Dec 5, 2024

simone-silvestri commented Dec 5, 2024

glwagner Dec 5, 2024

glwagner Dec 5, 2024

simone-silvestri Dec 5, 2024

glwagner Dec 5, 2024

glwagner Dec 5, 2024

glwagner Dec 5, 2024

simone-silvestri Dec 5, 2024

glwagner Dec 5, 2024

glwagner Dec 5, 2024

glwagner Dec 5, 2024

glwagner Dec 5, 2024

glwagner Dec 5, 2024

simone-silvestri Dec 5, 2024

glwagner Dec 5, 2024

simone-silvestri Dec 9, 2024

glwagner commented Dec 5, 2024

simone-silvestri commented Dec 9, 2024 •

edited

Loading

glwagner commented Dec 9, 2024

	q₁ = AtmosphericThermodynamics.vapor_specific_humidity(ℂ, 𝒬₁)
	q₀ = AtmosphericThermodynamics.vapor_specific_humidity(ℂ, 𝒬₀)
	Δq = q₁ - q₀

A first stab at a diagnostic surface temperature implementation #278

A first stab at a diagnostic surface temperature implementation #278

Conversation

simone-silvestri commented Nov 26, 2024 • edited Loading

simone-silvestri commented Nov 26, 2024 • edited Loading

simone-silvestri commented Nov 26, 2024 • edited Loading

simone-silvestri commented Nov 26, 2024

simone-silvestri commented Nov 26, 2024 • edited Loading

codecov bot commented Dec 2, 2024 • edited Loading

Codecov Report

simone-silvestri commented Dec 5, 2024

glwagner commented Dec 5, 2024

glwagner commented Dec 5, 2024

simone-silvestri commented Dec 5, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

glwagner commented Dec 5, 2024

simone-silvestri commented Dec 9, 2024 • edited Loading

glwagner commented Dec 9, 2024

simone-silvestri commented Nov 26, 2024 •

edited

Loading

simone-silvestri commented Nov 26, 2024 •

edited

Loading

simone-silvestri commented Nov 26, 2024 •

edited

Loading

simone-silvestri commented Nov 26, 2024 •

edited

Loading

codecov bot commented Dec 2, 2024 •

edited

Loading

simone-silvestri commented Dec 9, 2024 •

edited

Loading